Transport driven by biharmonic forces: impact of correlated thermal noise 
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We study an inertial Brownian particle moving in a symmetric periodic substrate, driven by a zero- 
mean biharmonic force and correlated thermal noise. The Brownian motion is described in terms of 
a Generalized Langevin Equation with an exponentially correlated Gaussian noise term, obeying the 
fluctuation-dissipation theorem. We analyse impact of non-zero correlation time of thermal noise 
on transport properties of the Brownian particle. We identify regimes where the increase of the 
correlation time intensifies long-time transport of the Brownian particle. The opposite effect is also 
found: longer correlation time reduces the stationary velocity of the particle. The correlation time 
induced multiple current reversal is detected. We reveal that thermal noise of non-zero correlation 
time can radically enhance long-time velocity of the Brownian particle in regimes where in the 
white noise limit the velocity is extremely small. All transport properties can be tested in the setup 
consisting of a resistively and capacitively shunted Josephson junction device. 
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I. INTRODUCTION 

Directed transport of particles can be driven by zero- 
average forces. The second principle of thermodynamics 
excludes it for systems at thermodynamic equilibrium. 
Therefore the system has to be at a nonequilibrium state 
and, moreover, some symmetries of the system have to be 
broken as e. g. the reflection symmetry of the spatially 
periodic potential or time- symmetry of the external force. 
Presence of both nonequilibrium conditions and asymme- 
try usually leads to transport and these two constituents 
form the ratchet concept. Transport in ratchet systems 
can be powered by mechanical, electrical, optical, chem- 
ical or electronic means. A large variety of models have 
been proposed and realized experimentally as solid state 
devices, cold atoms in optical lattices, superconducting 
devices, geometrically asymmetric lattices, to mention 
but a few, see the review pQ. 

The zero-average force can be deterministic or stochas- 
tic. An example of such a deterministic force is a time- 
periodic AC driving consisting of one or several har- 
monics. In some cases there is possibility of generating 
a charge or particle DC current from pure AC-driving, 
without the presence of an explicit static bias. The well- 
known phenomenon is harmonic mixing [2]. A paradig- 
matic model of this phenomenon is based on the Newton- 
Langevin equation which describes a particle moving in 
a spatially periodic potential (which mimics e.g. the mo- 
tion of a charge particle in a crystal) and driven by bi- 
harmonic force F(t) of two components of frequencies uj\ 
and 002 and amplitudes F\ and F2, respectively, 

F(t) = F 1 sin(a;it) + F 2 sm(uj 2 t + 0), (1) 

where <p is the phase-lag of two signals. Transport in such 
systems has been studied in various contexts, mainly in 
the overdamped regime [3], for moderate damping [4], 
also in dissipative quantum systems [5 , both experimen- 
tally and theoretically for cold atoms in the optical lat- 
tices [6 - 8 , and for driven Josephson junctions [9]. 



The sketch of the paper is as follows. In Sec. II 
we briefly present main elements of modelling for typ- 
ical systems exhibiting the ratchet effects when thermal 
fluctuations are described by Gaussian white noise. Sec. 
Ill contains the formulation of the corresponding model 
of non-Markovian dynamics for systems driven by bihar- 
monic signals and correlated thermal noise. In Sec. IV we 
analyse in detail influence of correlation time of thermal 
noise on transport properties. Finally, Sec. V provides 
summary. 

II. PARADIGMATIC MODEL 

An archetype of a ratchet is based on the Newton- 
Langevin equation 

mx + jx = -V'(x) + G(t) + f (t), (2) 

where m is a mass of the particle, 7 is the friction coeffi- 
cient and V(x) is a spatially periodic potential of period 
L, i.e. V(x) = V(x + L). The deterministic force G(t) 
is a non-biased function of zero average over some time- 
interval. The last term £(t) is a stochastic force which 
can model thermal equilibrium fluctuations and/or non- 
equilibrium noise. 

Variations of this equation are unlimited. Putting for- 
mally m = 0, the overdamped system is modelled. In the 
case m/0, inertial effects can be investigated. The po- 
tential V(x) can be symmetric or asymmetric, the force 
F(t) can be symmetric or not, the same for unless it 
mimics thermal noise, which is always symmetric. Prob- 
lem of symmetry of the stochastic driving £(£) is discussed 
inRef. [10]. 

If the process £(t) describes thermal equilibrium fluctu- 
ations then £(t) is zero- mean, Gaussian white noise with 
the Dirac delta auto-correlation function 

{Z(t)Z(s)) = *yk B To5(t-s), (3) 

where ks the Boltzmann constant and To denotes the 
temperature. Note that this process in not correlated 
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and in consequence its correlation time r c = 0. The noise 
intensity factor 27/^X0 follows from the fluctuation- 
dissipation theorem [11] . Gaussian white noise generates 
mathematically tractable models. In many real systems 
thermal noise is approximately white, meaning that the 
power spectral density is nearly equal throughout the fre- 
quency spectrum. Additionally, the amplitude of the sig- 
nal has very nearly a Gaussian probability density func- 
tion. However, an infinite-bandwidth white noise signal 
is purely a theoretical construction. By having power at 
all frequencies, the total power of such a signal is infi- 
nite and therefore impossible to generate. In practice, 
however, a signal can be "white" with a flat spectrum 
over a defined frequency band and it is a good approxi- 
mation of many real-world situations. E.g. in electrical 
systems, the white noise modelling is correct at any prac- 
tical radio frequency in use (i.e. frequencies below about 
80 GHz). In the most general case, which includes up 
to optical frequencies, the power spectral density of the 
voltage across the resistor depends on frequency and the 
white noise approximation fails. 

In the paper, we study transport in a spatially symmet- 
ric substrate and the corresponding potential is reflection 
symmetric, i.e. V(xq + x) = V(xq — x) for some fixed 
value xq. One of the simplest form of such a potential is 
given by the sinusoidal form 

V(x) = AVsin(27rx/L) (4) 

of the period L and the barrier height 2AV. The non- 
equilibrium driving G(t) is a biharmonic force and is cho- 
sen as a particular case of Eq. ([TJ , namely, 

G(t) = A[sm(Slt) + e sm(2Qt + (/))] . (5) 

Here A is the amplitude of the first harmonic, the factor 
e scales the second harmonics, so that it has the resulting 
amplitude eA. The angular frequency ft determines the 
time period T = 2tt/CI of the driving and <j) controls the 
phase shift between two components of the biharmonic 
signal ([5]), see Fig. [ljfor details. 

The Langevin equation ([2| with the potential Q has 
a similar form as an equation describing dynamics of the 
phase difference = ty(t) between the macroscopic wave 
functions of the Cooper pairs on both sides of the re- 
sistively and capacitively shunted Josephson junction, 
which is well known in the literature as the Stewart- 
McCumber model [12] [13] , namely, 

\2eJ \2e) R 2e 

= to m + *V-t m - (6) 

The left hand side is the total current through the junc- 
tion. It is a sum of three additive current contributions: a 
displacement current accompanied with a capacitance C 
of the junction, a normal (Ohmic) current characterized 
by the normal state resistance R and a Cooper pair tun- 
nel supercurrent characterized by the critical current Jo- 



in the right hand side, I(t) is an external current. The 
Johnson noise plays a role of thermal equilibrium 
noise associated with the resistance R and is modelled 
by Gaussian white noise of zero mean and the intensity 
taken according to the fluctuation-dissipation theorem. 
It is an evident correspondence between two models: the 
coordinate x = \£ — 7r/2, the mass m = (h/2e) 2 C the 
friction coefficient 7 = (h/2e) 2 (l/ R), the barrier height 
is determined by the relation AV = {h/2e)I$ and the 
period L = 2tt. The biharmonic signal G(t) in Eq. ^ 
corresponds to the external current I(t) . The velocity 
v = x corresponds to the voltage V across the junction. 
This correspondence allows for testing all transport prop- 
erties of the system ([2| and its generalizations by the 
setup consisting of a resistively and capacitively shunted 
Josephson junction device. 

III. NON-MARKOVIAN DYNAMICS 

If the correlation time r c of thermal noise is much 
smaller than the smallest characteristic (deterministic) 
time r x in the considered system then the white noise 
modelling is a good approximation. However, there are 
systems where this relation is not satisfied: the noise cor- 
relation time can be of order or greater than r x [T4HT6] . 
In this case we shell not assume the Markovian dynam- 
ics but rather invoke the Generalized Langevin Equation 
(GLE) for the time evolution of the system 

When the correlations comes into play the dynam- 
ics becomes non-Markovian and exhibits the memory- 
friction described by the GLE in the form [TTJ EG3 El CLE] 

mx{t)+ [ K(t- s)x(s) ds = -U'(x(i),t) + (7) 
Jo 

where £(£) is a stationary, zero-mean Gaussian stochastic 
process modelling correlated thermal noise and the full 
potential takes the form 

U(x, t) = V(x) - A[sm(nt) + e sin(2ftt + </>)]x. (8) 

The auto-correlation function C(t — s) of the correlated 
Gaussian thermal noise is related to the memory 
kernel K(t) via the fluctuation-dissipation relation [TTJ 

nana eei : 

C(t -s) = = k B T K(\t - s\). (9) 

Let us note that for more realistic modelling based on 
correlated thermal fluctuations we have to pay the price: 
Eq. Q is an integro-differential equation with a Gaus- 
sian random term which is difficult to handle by an- 
alytical or numerical means. The resulting process is 
neither Gaussian nor Markovian. In a general case, we 
do not know even a master equation for the single-event 
non-Markovian probability p(x, x, t) of the process deter- 
mined by Eq. ([7]). However, the problem is much easier 
to treat if the non-Markovian process x(t) is a projec- 
tion of a higher-dimensional Markovian process. In other 
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FIG. 1: (color online) Biharmonic external ac driving C(t) = 
A [sin (fit) + £sin(2fi£ + 0)] plotted for amplitude A — 1 and 
frequency fi = 1, for four different second harmonic amplitude 
e = 0.1,0.5,1,2 and four phase shifts = 0, tt/2, 1.65, 2.41 
which corresponds to the later commented figures (see [3] and 
[I] for details). 



words, it can be embedded in a higher (but finite) dimen- 
sional phase space and then Eq. ([7| can be converted to a 
set of first-order differential equations of Langevin type. 
It can be done if the kernel K(t) obeys an ordinary dif- 
ferential equation with constant coefficients. 



A. Correlated thermal noise 

Gaussian thermal noise is completely defined by 
its correlation function C(t). White noise corresponds to 



([3]). Examples of correlated noise are: the exponentially 
correlated noise [T9] , harmonic noise [20] , algebraically 
correlated noise [2TJ [22]. Other variations of noise are 
also considered, as e.g. when C(t) is a sum of exponen- 
tials [221124]. 

The simplest correlated noise is exponentially corre- 
lated noise (because of the smallest number of parame- 
ters) known as the Ornstein-Uhlenbeck (O-U) stationary 
stochastic Markov process for £(£) [15]. Its correlation 
function reads 

jk B T 



(fflt(8))=k B T K(\t-8\) 



-\ts\/7 



(10) 



Here r c is the correlation time of the O-U process. When 
the correlation time r c — » the correlation function ( [10] ) 
tends to the correlation function (|3l and the O-U process 
tends to white noise. Then Eq. (jT]) reduces to Eq. 
Below, we study the case of O-U noise. 



B. Markovian embedding dynamics 

For the case of the O-U noise, the integral kernel 
([9| is an exponential function of the time. The expo- 
nential function obeys a first-order differential equation 
with a constant coefficient. Therefore we can convert 
the integro-differential equation ^ into a set of ordinary 
stochastic differential equations. To this aim we define 
the auxiliary stochastic process 



w(t) 



1 

Tr 



-^- S )/ Tc i( 5 ) ds, 



(11) 



which is an integral part of Eq. Q. By means of the 
above relation the GLE ^ can be rewritten in the form 



mv{t) 
x(t) 

w(t) 



-U'(x(t),t) 
v(t), 

—w(t) + - 



v(t), 



£(*) 



-^2 7 k B T r(t), 



(12) 
(13) 

(14) 
(15) 



The last equation ( [15] ) corresponds to the O-U process 
with the exponential correlation function in Eq. ( 10 ) [25 . 



The stationary stochastic process T(t) describes Gaus- 
sian white noise of zero mean and correlation function 
(r(t)r(s)) = S(t — s). So, we embedded a non-Markovian 
process in a 4-dimensional space in which the process is 
Markovian. It is not a minimal dimension. It can be fur- 
ther reduced to the 3-dimensional space. To do it. let us 
note that for the linear combination z(t) = £(t)—w(t) we 
are able to subtract the two last relations and reproduce 
the three coupled Langevin equations [26 , namely, 



x(t) 
v(t) 

m 



v(t), 

1 



1 



m m 

Tr Tr T 



(16) 
(17) 

-y/2jk B T T(t). (18) 
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In the following we shall analyse this set of three coupled 
equations. 



C. Dimensionless variables 

The natural length scale for the system is settled by 
the period L of the potential V(x). For the adequate 
time scale we have to consider several time scales [27] . 
Here we define the characteristic time as Tq = mL 2 / AV 
which is convenient for studying inert ial effects. It can 
be obtained from Eq. Q by comparing the inertial term 
(with mass m) to the potential force V'{x) and inserting 
in both sides the characteristic quantities, for detail see 
Ref. [27 j. It leads to the scaling used throughout this 
paper 



X 



t 



t 



V T 

and finally to the rescaled evolution equations 



(19) 



X 
Y 



Y, (20) 
-W'{X) + a[sm(wt) + e sin(2wi + <p)} + Z, (21) 



1 



(22) 



where the explicit form of remaining coordinates reads 

L 



Y=-v, 



Z =AV Z - 



(23) 



The dot and prime denote the differentiation with re- 
spect to the scaled time t and the argument respectively. 
Rescaled velocity is designated by Y(t) and Z(t) sym- 
bolizes the dimensionless random force. The remaining 
re-scaled parameters are: 

1. the friction coefficient 7 = (7/771)70 = tq/tl is de- 
fined by the ratio the two characteristic times of 
the GLE - previously defined To and the relaxation 
time of the velocity degree of freedom tl =771/7, 

2. the correlation time f c = t c /to, 

3. the potential W(X) = V(x)/AV = W(X + 1) = 
sin(27rX) possesses the unit period and barrier 
height AV = 2, 

4. the amplitudes of the external deterministic stim- 
ulus scales as a = LA/ AV and the frequency 
uj = Qto (or the equivalent period T = 2tt/oj), 

5. the zero-mean white noise £(£) is correlated as 
(£(t)£(s)} = S(t — s) with a re-scaled noise intensity 
D = ksTo/AV and can be interpreted as the ratio 
of two energies, thermal energy and half of barrier 
height of the potential V(x). 

From now on we will use only above dimensionless vari- 
ables and shall omit the "hats" in all quantities of Eqs. 
d20b-d22). 



D. Method of analysis 

The analytical methods to handle nonlinear Brown- 
ian equations with memory friction are unknown to our 
knowledge. Therefore we will explore the peculiarities of 
the system (20)-(22) by means of numerically calculated 



long-time transport characteristics. In particular we shall 
focus on the current defined by the long-time averaged 
velocity v = (Y). The averaging is performed in the fol- 
lowing way: First the numerical mean over at least 10 3 
realizations of the GLE is calculated. This yields a time 
dependent quantity which we next average over one pe- 
riod of the external field T. In the process of calculating 
averages we must make sure that we initiate the process 
evolutions with unbiased initial conditions since the sim- 
ulated asymptotic, long time dynamics is not necessarily 
ergodic. We choose all initial positions and velocities to 
be uniformly distributed over one potential period [0, 1] 
and the interval v G [—2, 2], respectively. From the tech- 
nical point of view we have employed Stochastic Runge- 
Kutta algorithm of the 2 nd order [28] with the time step 
of [10 -3 -T- 10 _4 ]T. All numerical calculations have been 
performed using CUDA environment on desktop comput- 
ing processor NVIDIA Tesla C1060 [29]. 



IV. IMPACT OF CORRELATED THERMAL 
NOISE ON TRANSPORT 

In the long time limit, the averaged velocity can be 
presented in the form of a series of all possible harmonics, 
namely, 



lim (X) = v + v u (i) + v 2u (i) 



t— >OG 



(24) 



where v is a dc (time-independent) component and v nuj (t) 
are time-periodic functions which time average over a ba- 
sic period are zero. In figures we present only the dc 
component v of the long time, averaged velocity. From 
the symmetry property of Eq. ^ it follows when v van- 
ishes, for details see [8, 30-33 . In particular, in Ref. [31] 
it has been shown that if the force F(t) in Eq. ([T]) has 
the form 



F(t) = €1 cos(qujt + + 62 cos(pcjt + (^2) 



(25) 



and (p, q) are two coprime integers such that p + q is 
odd, the asymptotic velocity can be approximated by the 
expression 



v = Bc^Cr, cos(p</>i — q(j)\ + 0q) 



(26) 



provided the amplitudes e\ and 62 are sufficiently small. 
The quantities B and 0q depend on the parameters of the 
model and uj but neither on the amplitudes nor on the 
phases. 

We remind that in the Hamiltonian regime Eq. Q 
reduces to the form 



mx(t) = -U'(x(t),t) 



(27) 
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FIG. 2: (color online) Influence of the second harmonic of the 
external force G(t) on transport properties of the system |7|) 
is presented. Dependence of the drift velocity on both the 
relative phase <fi (horizontal axis) and the correlation time r c 
of thermal fluctuations (vertical axis) is depicted for increas- 
ing values of the amplitude e of the biharmonic driving G(t) 
(top to bottom). The other rescaled parameters read a = 4.2, 
u 4.9, 7 0.9, d 0.001. 



and for small amplitudes the dc velocity takes the sine- 
like form H |34] 



v oc sin ( 



(28) 



Mechanism of the transport generation in Hamiltonian 
systems is explained in Ref. [35]. The Hamiltonian 
regime can be realized in two case: (i) in the dissipa- 
tionless regime when 7 = and D = 0, (ii) in the limit 
of long correlation time, r c >> 1. 

For white noise case and in the weak damping regime, 
the dissipation-induced phase lag 0o occurs, i.e. [4[ [34] 



v oc sin 



{<t> ~ 0o), 



(29) 



where the phase lag 0o is determined by dissipation and 
vanishes in the Hamiltonian limit. 

We first address the issue of whether, and to which ex- 
tent, the non-zero correlation time r c of thermal fluctu- 
ations can influence transport properties. We thus start 
our analysis by studying the asymptotic dc velocity v 
as a function of both the time-symmetry breaking phase 
<j) and the correlation time r c . Our results are shown 
in Fig. 2 for selected values of the driving amplitude e 
of the second harmonics. A general conclusion from all 
cases presented in this figure is the destructive influence 
of the strongly correlated noise: for long correlation time 
r c the dc velocity is much smaller (virtually zero) than 
for the small-to-moderate correlation time. Let us note 
that the limit of the long correlation time (r c >> 1), cor- 
responds to the Hamiltonian regime in which transport is 
non-effective in comparison with the dissipative regime. 
A more accurate inspection reveals regions of weakly cor- 
related noise where reach diversity of transport charac- 
teristics can be observed. To identify them we search the 
section = const, of a constant phase in Fig. 2. The 
results are depicted in Fig. 3 where the average velocity 
is plotted as a function of the correlation time for three 
selected values of the amplitude e = 0.1,1,2 (the same 
as in panels a, c and d of Fig. [2]). The phase is fixed 
at the values <j) = 1.65 and (j) — 2.41. In Fig. 3, one 
can identify the so termed current reversal phenomenon: 
the velocity changes its sign as one of the parameters is 
varied. Here, even the multiple current reversal can be 
observed by changing the noise correlation time. E.g. 
for the amplitudes e = 1 in Fig. 3(a), the average ve- 
locity is negative for white noise, i.e. when r c = 0. If 
the correlation time r c increases, the velocity approaches 
zero value: the particle does not exhibit directed trans- 
port. Upon further increasing of r c , the velocity starts 
to increase to a positive- valued local maximum. Next, it 
again starts to decrease reaching a negative- valued local 
minimum. For further increase of r c , the velocity tends 
again to a positive-valued local maximum. Finally, with 
increasing r c it monotonically decreases towards smaller 
and smaller positive values. 

As the next point of analysis, we ask about the positive 
role of weakly correlated noise, i. e. whether small r c can 
intensify transport by enhancing the absolute value of 
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FIG. 3: (color online) The stationary averaged velocity vs 
the correlation time r c of thermal noise for three values of the 
relative amplitude e of the second harmonics and two values of 
the relative phase <j). Other parameters are: a — 4.2, uj — 4.9, 
7 0.9, d 0.001. 



the average velocity. The answer is yes. As an example, 
consider the case e = 2 in Fig. 3(b). For r c = 0, the 
velocity is negative and its absolute value is relatively 
large. If r c starts to increase, the absolute value of the 
velocity also increases attaining a global maximum at 
some fixed value of r c (of order 10 _1 ). In this example, 
the weakly correlated noise enhances transport in some 
interval of r c . On contrary, for the same amplitude e = 2 
but different value of the phase <p = 1.65, if r c starts 
to increase from zero, the absolute value of the velocity 
decreases, see the case e = 2 in Fig. 3(a). 

The influence of temperature is depicted in Fig. |4j 
where two various regimes are illustrated. In panel (a) 
we show the same curve as in panel (a) of Fig. 3 (the 
case e = l and do = 0.001). In this regime we can 
note the most pronounced influence of temperature in 
vicinity of the first maximum: if temperature increases, 
the first maximum is significantly reduced and velocity 
even changes its sign. So, by changing the temperature 
one can obtain the current reversal, see the curves for 
do = 0.001 and do = 0.01 in panel (a). We have searched 
a wide part of the parameter space in order to locate re- 
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FIG. 4: (color online) The stationary averaged velocity vs the 
correlation time r c of thermal noise for three values of the 
relative temperature do = 0.0001,0.001,0.01 and two sets of 
pairs e = 1,0= 1-65 (panel a) and s = 0.5, (j) = 3.24 (panel b). 
Please note the destructive role of the small correlation time 
r c on panel (a) and the antagonistic effect of the constructive 
impact of r c in panel (b). For large values of r c the average 
velocity is always much smaller and virtually is zero. Other 
parameters are: a — 4.2, s = 1, u = 4.9, 7 — 0.9. 



gions of temperature robustness. The conclusion is that 
for large values of the correlation time (in most cases 
r c > 1) the system is resistant to increase of temper- 
ature, at least to values presented in Fig. [4j On the 
contrary, for small correlation times, there are regimes 
where the system is sensitive to changes of temperature 
and increase of temperature responses in smaller values 
of the absolute value of velocity, see panel (b) of Fig. 
[4] It is a region of destructive role of temperature but 
the constructive role of short correlation times of ther- 
mal fluctuations: a first increase of r c from zero results 
in increase of the velocity. For moderate and long corre- 
lation time, the temperature smooths the dependence on 
the correlation time similarly as in both panels of Fig. [4] 
Finally, in Fig. [5] we present a regime where the am- 
plitude of the second harmonics is much greater than the 
first harmonics (here 10 times greater). In such a case, 
the degree of symmetry breaking is very small: the in- 
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FIG. 5: (color online) The stationary averaged velocity vs 
the correlation time r c of thermal noise for three values of 
the relative phase <j) — tt/2 blue (solid) line, tt green (dashed) 
line, 37r/2 red (dashed-dotted) line. Other parameters are: a 
0.4, s = 10, uj = 4.9, 7 0.9, d 0.001. 



second one (remember that if any of two harmonics is ab- 
sent, there is no transport in the system). Therefore the 
dc velocity is expected to be extremely small (practically, 
it is zero). It is true but only for short correlation times 
(from to 1 in Fig. 5). When we, however, increase r c a 
little bit more and pass the value of about 1, the signifi- 
cant net transport occurs. For some values of the phase 
(j) the velocity is negative velocity ((j) = tt/2); for other 
values of (j) it shows the positive response (0 = 3tt/2). 
In some other cases it exhibits the velocity reversal phe- 
nomenon as e.g. for <j) = tt. 

We have mainly concentrated the analysis on impact of 
the correlation time of thermal noise on transport prop- 
erties of the Brownian particle. Other parameters modify 
directed movement as well. We do not present all vari- 
eties but yet we refer to our webpage [36] where the more 
detailed analysis is presented. 



V. SUMMARY 




FIG. 6: (color online) Illustration of role of the thermal noise 
correlations and the relative phase of driving two harmonics 
on transport of the Brownian particle. It is a regime where 
long correlations of thermal noise destruct transport. The 
parameters: a — 4.2, e — 0.5, uj — 4.9, 7 — 0.9, d — 0.001 



fluence of the first harmonics is much smaller than the 



The richness and diversity of influence of non-zero cor- 
relation time of thermal noise on non-Markovian dynam- 
ics determined by Eq. ([7]) is summarized and visualized 
in Fig. 6. This picture looks like a landscape after in- 
tensive eruption of volcanoes: there are mountains with 
sharp summits and volcanoes with deep craters. This 
landscape can continuously be deformed to other compli- 
cated landscapes by changing parameters of the model. If 
you are long-suffering, you can discover your own original 
landscapes of non-Markovian dynamics driven by bihar- 
monic signals. The diversity of subtle structures hidden 
in Eq. ^ is infinitely large. It comes to our mind a 
loose analogy to a Mandelbrot set which conceals beau- 
tiful forms of Julia sets, fractals and extremely difficult 
mathematics. 
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